%  5 BY, nT (GSM)
%  6 BZ, nT (GSM)
%  7 Flow pressure, nPa
%  8 Electric field, mV/m 
% 
load c:\manoj\projects\substorm\omni_2000_2005_1_min data;
L =data(:,5) == 9999.99;
data(L,5) = NaN;
L =data(:,6) == 9999.99;
data(L,6) = NaN;
L =data(:,7) == 99.99;
data(L,7) = NaN;
L =data(:,8) == 999.99;
data(L,8) = NaN;
fday_ace = datenum(data(:,1),0,data(:,2),data(:,3),data(:,4),0)- datenum(2000,1,1) + 18/(24*60);
% 18 minutes delay go the maximum picking 
load c:\manoj\projects\ace\substorm;

load c:\manoj\projects\substorm\substorm_omni_index substorm_omni_index substorm_omni_centerpoint;

component = 6;% 5- by, 6- bz



% for i = 1:100,
%     
%     L = fday_omni_1_min_2000 > fday_substorm(i)-radius_epoch & ...
%                fday_omni_1_min_2000 < fday_substorm(i)+radius_epoch;
%     
%     if sum(L) > 1,
%     seltime = fday_omni_1_min_2000(L);nci
%     seldata = data(L,:);
%     
%     plot(seltime,seldata(:,5));
%     hold on;
%     plot(seltime,seldata(:,6));
%     plot(fday_substorm(i),0,'r*');
%     hold off;
%     pause;
%     end;
% end;

%%
radius_epoch = 40/(24*60);

n_counter = 0;
n_increasing = 0;
n_buildup = 0;
n_Nturning = 0;
n_sustained = 0;
n_other = 0;
n_bznegative = 0;
n_bzflat = 0;


flag_counter = 0;
flag_buildup = 0;
flag_Nturning = 0;
flag_sustained = 0;
flag_increasing = 0;
flag_bznegative = 0;
flag_bzflat = 0;




%%
clear selected_substorm;
for i = 1:length(fday_substorm),
    %find the nearest ACE data for the substorm
%     [trash ii] =  min(abs(fday_ace - fday_substorm(i)));
    
       
    %test whether it is near enough
    
%     if trash < 1/(24*60),
        
        %test at least 70 % of data are available within +/- 40 minutes
        
%             L =fday_ace >= fday_substorm(i)-radius_epoch & ...
%                fday_ace <= fday_substorm(i)+radius_epoch;
           
         %use saved indices
         
           selected_data = data(substorm_omni_index(i,1):substorm_omni_index(i,2) , component);
           selected_time = fday_ace(substorm_omni_index(i,1):substorm_omni_index(i,2));
           

           %saving indexes
           
%            index = find(L==1);
%            
%            substorm_omni_index(i,:) = index([1,end]);
           center_point = find(selected_time == fday_ace(substorm_omni_centerpoint(i)));
%            substorm_omni_centerpoint(i) = ii;
           
           if sum(~isnan(selected_data)) >= floor (0.7 * radius_epoch * 2 * 24 * 60),
               
                   % test the imf -ve condition prior to substorm onset
                   
                   if sum(selected_data(center_point-29:center_point) < 0) >= 22,
                     n_buildup = n_buildup + 1;
                     flag_buildup = 1;
                   end;
                    
                       % test rapid turning of imf bz
                       
                   if diff(selected_data(center_point:center_point+1)) > 0.375,
                      n_Nturning = n_Nturning + 1;
                      flag_Nturning = 1;
                   end;
                   
                           %test sustained turining
                   b = polyfit(selected_time(center_point+1:center_point+10),selected_data(center_point+1:center_point+10),1);
                   if b(1)/(24 * 60)  >= 0.175,
                       n_sustained = n_sustained + 1;
                       flag_sustained = 1;
                   end;
                      
                       
                  %fprintf('i = %d , fday substorm = %6.2f slope = %5.4f\n', i, fday_substorm(i), b(1)/(24 * 60));
                  
                  if (all(selected_data(center_point+1:center_point+3) >= selected_data(center_point) + 0.175) && ...
                      all(selected_data(center_point+4:center_point+10) >= selected_data(center_point) + 0.45)),
                                n_increasing = n_increasing + 1;
                                flag_increasing = 1;
                  end;
                  
                  %if imf bz was negative for most of the time around onset
                if sum(selected_data(center_point-29:center_point+29) < 0) >= 45
                    n_bznegative = n_bznegative + 1;
                    flag_bznegative = 1;
                end;
                   
                
                %find if IMF bz was almost flat and have low scatter
             
                   [rs1,rstat] = robustfit(selected_time,selected_data);
                 %  if abs(rs1(2)/(24 * 60))  < 0.08  & rstat.ols_s < 0.6,
                   if abs(rs1(2)/(24 * 60))  < 0.03  & rstat.ols_s < 0.4,
                       n_bzflat = n_bzflat + 1;
                       flag_bzflat = 1;
                   end;

                 
                % if (flag_buildup & flag_Nturning & flag_sustained & flag_increasing),
                if (flag_bzflat),
                      
                      n_counter = n_counter + 1;         
                      selected_substorm(n_counter) = fday_substorm(i);
                      fprintf('i = %d\n',i);        
                else
                    
                    n_other = n_other + 1;
                    non_selected_substorm(n_other) = fday_substorm(i);
                         
               end;
                           
                           
                            flag_increasing = 0;
                            flag_buildup = 0;
                            flag_Nturning = 0;
                            flag_sustained = 0;
                            flag_bznegative = 0;
                            flag_bzflat = 0;

                           
 
    end;
end;
       

%%                           
                       
n_counter = 0;

[dd,ii] = intersect(fday_substorm,selected_substorm);
%[dd,ii] = intersect(fday_substorm,non_selected_substorm);
for i = 1:length(ii),
k = ii(i);   
if size(substorm_omni_index(k,1):substorm_omni_index(k,2),2) == 80,
n_counter = n_counter + 1;
A(n_counter,:) = substorm_omni_index(k,1):substorm_omni_index(k,2);
end;
end;

plot(-radius_epoch*1440+0.5:radius_epoch*1440-0.5,((reshape(data(A,6),[size(A,1),size(A,2)])))');
hold on
plot(-radius_epoch*1440+0.5:radius_epoch*1440-0.5,nanmean((reshape(data(A,6),[size(A,1),size(A,2)]))),'g','LineWidth',4);
title(sprintf('Stack plot of IMF Bz nT N = %d for slope LE 0.03 nT/min AND sigma < 0.4 nT',n_bzflat));
ylabel(' IMF Bz nT');
xlabel('Time in minutes from substorm onset');
 

%%

% n = 0;
% for i = 1:length(fday_substorm),
% 
% selected_bz = data(substorm_omni_index(i,1):substorm_omni_index(i,2) , 6);
% selected_by = data(substorm_omni_index(i,1):substorm_omni_index(i,2) , 5);
% selected_pressure = data(substorm_omni_index(i,1):substorm_omni_index(i,2) , 7);
% selected_ey = data(substorm_omni_index(i,1):substorm_omni_index(i,2) , 8);
% selected_time = fday_ace(substorm_omni_index(i,1):substorm_omni_index(i,2));
% 
% if sum(isnan(selected_bz)) < 20,
%  n = n+ 1;
%  
% 
% [rs1,rstat] = robustfit(selected_time,selected_bz);
% 
% %  fit_data(n,1) =  rs1(2)/(24*60);
% %  fit_data(n,2) = rstat.ols_s;
% 
% if abs(rs1(2)/(24 * 60))  < 0.04 & rstat.ols_s < 0.8,
% subplot(311);
%     plot(selected_time,selected_bz,'r');
%     hold on;
%     plot(selected_time,selected_by,'g');
%     hold off;
%     
% title(sprintf('slop= %5.3f rstat = %6.3f', rs1(2)/(24*60), rstat.ols_s));
% axis([-inf inf -10 10]);
% 
% subplot(312);
% plot(selected_time,selected_pressure);
% axis([-inf inf -2 2]);
% 
% 
% subplot(313);
% plot(selected_time,selected_ey);
% axis([-inf inf -10 10]);
% 
% 
% pause;
% 
% end;
% 
% end;
% 
% end;
% 
% 
%         
% 
% 
%         
%     
%     
%     
%     